/*
 * 单原子自旋可以向上、向下，磁矩为mu，外场B向上，哈密顿量为
 *  H = - mu B Sz
 *  下面的代码输出马尔科夫链在屏幕上
 */

#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>

void markov_256atom(int n, double muB, double kT, double & ave_e, double & ave_sz){

	double sum_E = 0, sum_Sz = 0;
	int sz[256]; for(int i=0;i<256;i++) sz[i] = 1;
	double E, Sz;
	for(int i=0;i<n;i++){
		Sz=0; for(int j=0;j<256;j++) Sz += sz[j]; 
		sum_Sz += Sz;
		sum_E += - muB * Sz;
		for(int j=0;j<256;j++){
			if(sz[j] ==1){
				if( (double)rand()/RAND_MAX  < exp(-2*muB/kT) )
					sz[j] = -1;
			}
			else{
				sz[j] = 1;
			}
		}
	}

	ave_e = sum_E /n/256;
	ave_sz = sum_Sz /n/256;
}

int main(){

	int n=1000;
	double muB = 2;
	double kT, ave_e, ave_sz;
	ofstream fp("mc_markov_demon_sampling_256atom.txt");
	for(kT = 0; kT <50; kT += 0.1){
		markov_256atom(n, muB, kT, ave_e, ave_sz);
		fp<< kT<<"\t"<<ave_e<<"\t"<<ave_sz<<endl;
	}
	fp.close();
	return 0;
}
